
%% runParamVariation.m 
% 
% This code loads a mode from an XLS-07 file. 
% A single user-defined parameter is perturbed on a grid. 
% The model is then solved for each parameter. 
% A number of theoretical moments and plots the IRFS for the observables are displayed.  
% Use to understand the role of a single parameter in shaping the transmission of shocks and moments. 
%
% See runParamSingle.m for a code that does this the same for a
% a single parameter
%
% Alejandro Justiniano (c) Oct 30 2013
%
% Updated: Todd Messer 9/30/2014
%   Will print
%=========================================================================
clear all; 
close all; 
clc; 

cucd=cd; 

%% INPUTS
% Function Handle 
funcHandle=@modBMJ_newBaseDepU; 
addsol.GDPCWCons=0; 
addsol.measureEQ.wage.numSeries     = 2;
addsol.measureEQ.wage.nameSeries    = {'LEPRIVA_Core' 'LS_Core'};
addsol.measureEQ.price.numSeries    = 3;		
addsol.measureEQ.price.nameSeries   = {'JCXFEBM_LD100'	'JCMXFEM_LD100' 'PCUSLFEN_LD100SA'};

% Name and location of the XLS-07 filename that contains the mode 
xls.filename='gdrift_0-50 10 06 14'; 
xls.sheet='ModePrecision'; 
xls.path=fullfile(cucd,'BMJ_new','BaseDepU','RemoveCons_ShutdownGDrift_M2','end 2007.50 csminwel gdrift_0-50 10 06 14'); 

% Model will be solved for each value of *paramVariation* in increasing
% order, adding to the grid below the value loaded in the mode 
paramVary.name='STDgdrift'; 
paramVary.grid=[1 2]; 

% Number of horizons for the IRFS 
dim.horIRF=16; 
% Number of horizons for the forecast error decompositions 
dim.horDec=20; 

% Shock to Plot
paramShock.name='Government Drift';

%==========================================================================
%% Load XLS File with Mode 
cd(xls.path);
[parVec,parNames]=xlsread(xls.filename,xls.sheet);
cd(cucd); 

%==========================================================================    
%% Solve Model 
% Preliminary Solution as Test
[GG,RR,CONS,eu,SDX,ZZ,~,~,flags,~,stateNames,shockNames]=feval(funcHandle,parVec,[],addsol); 
if ~isequal(eu,[1;1])
    error('No Solution')
end

% Extract the names of the shocks
if isstruct(shockNames)==true 
    shockNames=shockNames.long; 
end 

% Extract the names of the observables as defined in ZZ 
[obsNames,~]=extractObsNames(ZZ(:,:,end),stateNames);  

% Sizes 
dim.y=size(ZZ,1);       % Observables 
dim.s=size(GG,1);       % Full state of transition equation 
dim.eta=size(SDX,1);    % Innovations 
dim.paramVary=numel(paramVary.grid)+1; % Variations  

%% Crate a matrix of parameters with the parameter perturbation 
paramVary.pos=cellposition(paramVary.name,parNames); 
paramVary.grid=sort( [parVec(paramVary.pos) paramVary.grid] ); 
paramVary.mat =repmat( parVec ,[1 dim.paramVary] ); 
paramVary.mat(paramVary.pos,:)=paramVary.grid; 
paramVary.cell=cell( dim.paramVary ,1 ); 
for ii=1:dim.paramVary 
    paramVary.cell(ii)={ sprintf( [char(paramVary.name),'=','%2.2f'],paramVary.grid(ii) )}; 
end 

%% Storage Matrices 
paramVary.stdObs   =nan(dim.y, dim.paramVary );     % Standard deviation of the Observables 
paramVary.stdStates=nan(dim.s, dim.paramVary );     % Standard deviation of the States 
paramVary.irfObs   =nan(dim.y, dim.horIRF+1, dim.eta , dim.paramVary ); % IRF of the Observables
paramVary.irfStates=nan(dim.s, dim.horIRF+1, dim.eta , dim.paramVary ); % IRF of the States
paramVary.EUOK=zeros(dim.paramVary); % Indicator if model has unique stable solution 

%% Evaluate Different Models
for varInd=1:dim.paramVary 
    
    [GG,RR,CONS,eu,SDX,ZZ]=feval(funcHandle,paramVary.mat(:,varInd),[],addsol); 
    if isequal(eu(:),[1;1]) 
        paramVary.EUOK(varInd)=1; 
    else 
        continue 
    end 
    
    [stdObs,stdStates,irfObs,irfStates,ACZeroObs,ACOneObs,VDecObs,VDecStates,...
        VDecHorObs,VDecHorStates]=...
        irfandthmom(GG,RR,SDX,ZZ,dim.horIRF,dim.horDec,zeros(size(ZZ,1),1),0);

    paramVary.stdObs(:,varInd)=stdObs; 
    paramVary.stdStates(:,varInd)=stdStates; 

    paramVary.irfObs(:,:,:,varInd)=irfObs; 
    paramVary.irfStates(:,:,:,varInd)=irfStates; 

end 

%==========================================================================
%% OUTPUT 
% 1. stdStates: std. dev. of the states (st)
printcell([crtcell(['Theoretical Std. Dev. States' paramVary.cell']); crtcell(paramVary.stdStates,stateNames(:))])
 
% 2. stdObs: std. dev. of the observables (z) 
printcell([crtcell(['Theoretical Std. Dev. Observables' paramVary.cell']); crtcell(paramVary.stdObs,obsNames(:))])

% 3. Asymptotic (inf horizon) decomposition of the observables 
printcell( crtcell( VDecObs',shockNames,obsNames,'Asymptotic Variance Decomposition')); 

% 4. Impulse Response Functions
fig.vec=zeros(dim.eta,1); 
fig.cols=3; 
fig.rows=ceil( dim.y / fig.cols ); 
fig.xAxis=(0:dim.horIRF); 
fig.xAxis=fig.xAxis(:); 
fig.colors=hsv(3);

% Plotting
shplot = cellposition(paramShock.name,shockNames);

for shind=shplot; %1:dim.eta
    
    fig.vec(shind)=figure('Visible', 'off','PaperPosition',[0.75 0.5 7.5 10]);
    
    for serind=1:dim.y 
        
        subplot( fig.rows, fig.cols, serind ); 

        for i = 1:dim.paramVary
            plot(fig.xAxis,squeeze( paramVary.irfObs(serind,:,shind,i)),'color',fig.colors(i,:),'LineWidth',1.5); 
            hold on
        end
        
        if serind==dim.y
            legend(paramVary.cell,'location','southoutside')
        end
        
        title(char(obsNames{serind})); 
        
        set(gca,'box','off'); 
        
        axis tight; 
        
    end

    suptitle(['Shock: ',char(shockNames{shind})],14);

end         

%% HARDCODED: Difference between ee_glogDrift - ylog
pos.g = cellposition('gov+nx drift (\epsilon^{gdrift})',stateNames);
pos.y = cellposition('GDP',stateNames);
diffShare = cell(3,1);
fig.vecShare=figure('Visible', 'on','PaperPosition', [0 0 11 8.5], 'PaperSize', [8.5 11],'PaperOrientation','landscape');
numVary = max(dim.paramVary,1);
for i = 1:numVary
    
    disp(5)
    
    diffShare{i} = squeeze(paramVary.irfStates(pos.g,:,shind,i)) - squeeze(paramVary.irfStates(pos.y,:,shind,i));
    
    subplot(2,2,1);
    plot(fig.xAxis,squeeze(paramVary.irfStates(pos.g,:,shind,i)),'color',fig.colors(i,:),'LineWidth',1.5);
    title('lnlin Government Drift')
    hold on
    
    subplot(2,2,2);
    plot(fig.xAxis,squeeze(paramVary.irfStates(pos.y,:,shind,i)),'color',fig.colors(i,:),'LineWidth',1.5);
    title('lnlin Output')
    hold on
    
    subplot(2,2,3);
    %     plot(fig.xAxis,diffShare{i} + log(0.17),'color',fig.colors(i,:),'LineWidth',1.5);
    plot(fig.xAxis,diffShare{i},'color',fig.colors(i,:),'LineWidth',1.5);
    title('lnlin Government Drift - lnlin GDP')
    hold on
    
    subplot(2,2,4);
    load lnDataShare
    plot(Dates(Dates>=1993 & Dates<2008),lnDataShare(Dates>=1993 & Dates<2008),'color','blue','LineWidth',1.5);
    title('Data: ln Residual (Y-G-I) - ln Output')
    hold on
end

subplot(2,2,3);
legend(paramVary.cell);
axis tight;



%% Save all as PDFs
cd(cucd)
foldername = ['OutputGraphs-' datestr(now)];
foldername = regexprep(foldername,':','-');
foldername = regexprep(foldername,' ','_');
mkdir(foldername)
cd(foldername)
for shind = shplot
    filename = regexprep(shockNames(shind),' ','_');
    print(fig.vec(shind),'-dpdf',[filename{1} '.pdf'])
    
end    
print(fig.vecShare,'-dpdf','GovernmentShare.pdf')
cd(cucd)
close all